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High-throughput ab-initio calculations, cluster expansion techniques and thermodynamic model¬ 
ing have been synergistically combined to characterize the binodal and the spinodal decompositions 
features in the pseudo-binary lead chalcogenides PbSe-PbTe, PbS-PbTe, and PbS-PbSe. While our 
results agree with the available experimental data, our consolute temperatures substantially improve 
with respect to previous computational modeling. The computed phase diagrams corroborate that 
the formation of spinodal nanostructures causes low thermal conductivities in these alloys. The 
presented approach, making a rational use of online quantum repositories, can be extended to study 
thermodynamical and kinetic properties of materials of technological interest. 

PACS numbers: 64.75.Qr, 71.15.Nc, 81.30.Bx 


I. INTRODUCTION 

For decades, the physical properties of lead chalco¬ 
genides have generated substantial interest in a number 
of fields, in particular for applications in semiconductor 
technology.^ PbS, PbSe, and PbTe have distinct struc¬ 
tural and electronic properties compared to III-V and II- 
VI compounds. These include high carrier mobilities, nar¬ 
row band gaps with negative pressure coefficients, high 
dielectric constants, and a positive temperature coeffi¬ 
cient. In addition, PbS, PbSe, and PbTe were predicted 
to be weak topological insulators, with a band inversion 
observed at the N point of the distorted body-centered 
tetragonal Brilllouin zone.® These important and varied 
properties have allowed lead chalcogenides to be used ex¬ 
tensively in optoelectronic devices such as lasers and de¬ 
tectors, thermophotovoltaic energy converters, and ther¬ 
moelectric materials. 

As thermoelectric materials, lead chalcogenides may 
exhibit electrical conductivities, a, in excess of 
2-4 • 10“^n“^cm“^, thermopowers, S, around 150 
/iVK“^, and thermal conductivities, k, on the order of 
1-2 Wm“^K“^. This leads to figures of merit ZT = 
o’S'^fK larger than 1 at high temperatures, T. Such 
outstanding performances are due on the details of the 
electronic structure,and on the ability to dramati¬ 
cally reduce the thermal conductivity with alloying and 
nanostructuring. While pure lead chalcogenides are 
attractive on their own, their alloys are even more in¬ 
teresting. The appeal arises from their mechanical and 


electronic tunability, which can be optimized for specific 
technological needs.For example, the PbTei_a;Sea; 
pseudo-binary system has a higher ZT value than its cor¬ 
responding binary forms.Thallium doping in PbTe 
causes changes in the electronic density of states, in¬ 
creasing the ZT value to 1.5 at 773 K}^ Similarly, a 
ZT of 1.3 at 850 K was reported for aluminum doped 
PbSe.^® Furthermore, Pb 9 . 6 Sbo. 2 Teio-a:Sba; is known to 
exhibit lower thermal conductivity and a higher ZT 
than PbSei_a;Tea;.^® This is also true for nanostructured 
(Pbo,95Sno.o5Te)o.92(PbS)o.o85 as the low thermal conduc¬ 
tivity leads to a ZT = 1.5 at 642 The properties 
of this group of pseudo-binaries depend greatly on the 
atomic details of the material’s morphology. This can 
be partially understood in terms of thermodynamical fea¬ 
tures. The excellent thermoelectric performances, for ex¬ 
ample, were tentatively ascribed to the limited miscibility 
of the components. This gives rise to structural inho¬ 
mogeneities that lower the thermal conductivity without 
damaging the electronic transport.^® Such control over the 
morphology could be also used to optimize functionalities 
associated with topological effects. 

In this work we study the phase diagram of lead chalco- 
genide pseudo-binaries. We predict quantitatively the 
boundary of the solid solution (the binodal curve that 
defines the region of miscibility), as well as the spinodal 
region. These features are key for rationalizing and hon¬ 
ing synthesis and characterization of optimized systems. 
To the best of our knowledge, our study is the first to 
completely and accurately report such characterization. 
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The phase diagram is essential for properly establishing 
manufacturing processes. There has been only one pre¬ 
vious attempt to model phase diagrams of lead chalco- 
genides using thermodynamic modeling (TM),^® which 
predicted consolute temperature {T^) values far from 
those reported in experimental studies.^® The disagree¬ 
ment was attributed to the difficulty of an exhaustive ex¬ 
ploration of the different configurations for each compo¬ 
sition. Here, we built upon the synergy between cluster 
expansion (CE) techniques, high-throughput (HT) ab ini¬ 
tio calculations,and thermodynamical modeling to 
find acceptable agreement between our Monte Carlo sim¬ 
ulations (MC) and the available experimental results. 


II. METHODOLOGY 


A. Thermodynamic Modeling 


Pseudo-binary systems are represented by the for¬ 
mula (AxjiBxg)aCc (or AxBi-xC) with mole fractions 
XA and xb of elements A and B respectively, related by 
xa+xb=^- The small letters a and c represent number 
of sites per formula.®®’®^ The Gibbs energy of such iso- 
structural pseudo-binary systems can be written as:®®’®^ 


G{A,B)aG^ = XaGa^Cc AXbGb^Cc 
-\-kBT{xAln{XA) + XBln{xB)) XaXbLa,B:C, 


( 1 ) 


where and Gb^^c^ represent the Gibbs free energy 

of AaCc and BqGc materials. These two variables can be 
computed at any temperature by fitting available exper¬ 
imental data®^ to the polynomial form®® shown in Equa¬ 
tion (2): 

G{T) = a + bT + cTlnT + dT^ + eT"^ + fT^ ■ (2) 

The third term in the Equation (1) is the entropy of mix¬ 
ing, and the last term is the excess energy of mixing that 
represents the non-ideality of the system. In contrast to 
the entropic term, the excess energy parameter can take 
negative or positive values. If La,B:C is negative, it in¬ 
dicates that the system tends to create a solid solution. 
A positive value of La.b-.c indicates a repulsive interac¬ 
tion between phases, penalizing formation of intermediate 
alloys. To find the excess mixing energy, the composition- 
dependent interaction parameter La,B:C can be calculated 
with Equation (3): 


AH = xaxbLa,b-.c- (3) 

The enthalpy of formation, AH, is defined as: 

ah = E(^a,b)^Cc ~ xaEa^Cc ~ xbEb^Cc^ (4) 

where E(^ab) CcJ Ea^Coj Eb^c^ sxe the total ener¬ 
gies of compounds (AB)a Gc, Aa Gc and Ba Cc, respectively. 
These energies can be found from the fully relaxed struc¬ 
tures using density functional theory (DFT). 


A combination of high throughput ab initio calcula¬ 
tions and thermodynamic modelling are used to predict 
the interaction parameter, The result is a zero 

temperature approximation of the actual value. The 
most common method to describe the composition depen¬ 
dent interaction parameter is the Redlich-Kister equation, 
33,34,39 |;]^g interaction parameter is written in a 

polynomial form: 

n 

i^O 

We fit this polynomial to the formation enthalpy data 
obtained from DFT calculations. We checked that an n = 
2 polynomial is enough to obtain a good fit to the data, 
so that only Lq, Li, and L 2 need to be determined: 

L(x) ~ Lo -I- xLi x^ L 2 (6) 

to compute the interaction parameter. 

The Gibbs free energy is composition and tempera¬ 
ture dependent. The main computational challenge lies 
in characterizing many configurations for many compo¬ 
sitions. Some authors have attempted to tackle the is¬ 
sue, by generating, few configurations and/or few com¬ 
positions. 28,36-38 jjg]-g^ have chosen to challenge the 
issue more drastically, by reling on the advantages of a hy¬ 
brid cluster expansion - high throughput approach^® fea¬ 
turing: exhaustive exploration of different configurations 
for different compositions (GE), minimization of compu¬ 
tational cost by reducing the number of the ab-initio cal¬ 
culations (GE), analysis of many compositions as typical 
of HT methods^®’"^®, and rational use of online repositories 
(AFLOWLIB.org)42. 

B. Cluster expansion 

In the cluster expansion technique, the configurational 
energy of an alloy as written as a sum of many-body oc¬ 
cupation variables {c}:^® 

E{<7) = Jq + Jjdj + -!-■■■, (7) 

i ij 

where Jq, Ji, Jij-, etc. are known as effective cluster inter¬ 
actions and must be determined. 

The above equation can be rewritten into symmetrically 
distinct sets of clusters, a: 

E{a) ='^maJa(W(Ji), (8) 

Oi i^OL 

where TOq, represents symmetrically equivalent clusters a 
in a given reference volume.The Jq, parameters are ob¬ 
tained by fitting a relatively small number of DFT calcu¬ 
lated energies. 

The reliability of the predicted energy may be deter¬ 
mined using the cross-validation score GV: 

1 ^ 

{GVf = -Y,{Es - Es)\ 

^ s=i 


( 9 ) 
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where Eg represents the calculated energies from DFT 
and Es is the predicted energies from CE. 

The enumeration of configurations, calculation of the 
effective interaction parameters, determination of ground 
state structures, and prediction of more structures was 
performed with the Alloy Theoretic Automated Toolkit 
(ATAT).^® Calculated phase diagrams were obtained with 
Monte Carlo (MC) simulations performed with phb code. 
The algorithm automatically follows a given phase bound¬ 
ary and is provided by the ATAT package. 


at low temperatures. At high temperatures close to Tc, it 
becomes convex, with one minimum. 

The binodal curve Tc{x) is defined by the horizontal 
tangent points of the Gibbs free energy, G. When T <Tc, 
the alloy starts decomposing. Additionally, our calcula¬ 
tions let us determine the spinodal curve that discrimi¬ 
nates metastable and unstable regions in the pseudo-alloy 
phase diagram. The spinodal curve is the locus of the 
points where the second derivative of G is equal to 0: 


C. High-throughput ab initio calculations 


d^G{AB)aC^ 

dx^ 


( 10 ) 


All DFT calculations were carried out by using the 
Automatic-Flow for Materials Discovery (AFLOW) ^2,48,49 
and DFT Vienna ab initio simulation program (VASP).^° 
Calculations were performed using AFLOW standards. 
We use the projector augmented wave (PAW) pseudopo- 
tentials^^ and the exchange and correlation functionals 
parametrized by the generalized gradient approximation 
proposed by Perdew-Burke-Ernzerhof.®^ All calculations 
use a high energy-cutoff, which is 40% larger than the 
maximum cutoff of all pseudopotentials used. Reciprocal 
space integration was performed using 8000 more k-points 
than the number of atoms. Spin-orbit coupling was not 
treated in the calculations due to its minimal influence in 
AH (smaller than 1.5 meV/atom). Structures were fully 
relaxed (cell volume and ionic positions) such that the en¬ 
ergy difference between two consecutive ionic steps was 
smaller than 10“^ eV. 

PbS, PbSe and PbTe crystallize in the NaCl structure 
and belong to the Fm3m space group (# 225). Space 
and point group symmetries of intermediate composition 
structures were determined using AFLOW. 

III. RESULTS AND DISCUSSION 
A. The PbSei_2:Te2, alloy 

In agreement with experimental data, we found that 
PbSe and PbTe are immiscible systems at 0 K. This is 
shown in Figure. 1(a), where formation enthalpies are pos¬ 
itive for all compositions (0 < a: < 1). The CE predicted 
energies are in excellent agreement with DFT calculated 
structures, with a cross validation score of 4-10“"^. Our 
quantitative results confirm the Hume-Rothery rules. 
These rules qualitatively predict the miscibility of two 
metals based on four properties: atomic radius, crys¬ 
tal lattice, valence and electronegativity. Amongst the 
chalcogens, the atomic radius changes from 1.04 A in S, 
to 1.17 A in Se, to 1.37 A in Te. These size variations cre¬ 
ate a mismatch between the lattice parameters of PbTe 
and PbSe, causing incoherence in the interface and phase 
decomposition, eventually. AG(x) diagrams at different 
temperatures are plotted in Figure. 1(b). It can be seen 
the A G function has two minima and a single maximum 


We computed first and second derivatives of the Gibbs 
free energy within our thermodynamical model. In order 
to obtain the Lq, Li, L 2 that are necessary to compute 
L{x) at any composition, we use highly symmetric struc¬ 
tures (HSS)®®’®® to fit Equation (3). HSS have a larger 
degeneracy and thus greater weight in the properties of 
the ensemble; particularly at high temperatures that are 
close to the spinodal decomposition. After computing the 
Ln constants, G(Se,Te)aPba i® obtained using Equation (1). 
The results obtained from the thermodynamic model are 
compared with our MC results, previous theoretical pre¬ 
dictions,^® and experimental data.^®’®^ 

Our calculations encompass the entire range of con¬ 
centrations, and reproduce the asymmetry observed ex¬ 
perimentally in the binodal curve. In order to quantita¬ 
tively compare all data, we analyzed the consolute tem¬ 
perature or upper critical solution temperature (Tc) de¬ 
scribing the lowest temperature at which both phases are 
miscible at any composition. Experimentally PbSei_a;Te 2 ; 
alloy presents, at a; = 0.4, a Tc closer to 623 This 

quantity is far from the value predicted by Boukhris et 
^ which is around 106 K. The consolute tempera¬ 
ture predicted by our thermodynamic model at 520 K 
(16.5% error) is around x = 0.34 (see Figure. 1(c)). Our 
results quantitatively improve the prediction of T^ with 
of Boukhris et al.^® Moreover, this value approaches the 
results obtained by much more expensive techniques such 
as MC, in which we obtain a value close to 550 K (11.7% 
error). Discrepancies between MC and TM at larger Te 
(Se) concentrations arise from difficulties in converging 
MC calculations in the dilute limit of Se (Te). The exper¬ 
imental miscibility gap presents a slight asymmetric form 
that is reproduced by MC and our TM. This asymmetry 
is observed in experiments, but was not seen in previous 
theoretical work.^®. This phenomenon will be discussed 
in the next section. 

In contrast to the binodal curve, the spinodal curve is 
quite symmetric. The combination of the symmetric spin¬ 
odal and asymmetric binodal curves at the Se-rich region 
causes nucleation at higher temperatures. This informa¬ 
tion is very important for fine tuning synthesis protocols 
to obtain the desired morphologies. 
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PbSe Xxe PbTe 

FIG. 1. (a) Formation enthalpies of the PbSe-PbTe structures 
using DFT calculations (O) CE technique (□). Highly 
symmetryic structures are represented by A and the fitting 
of these points to obtain the interaction parameter is plotted 
with a blue dashed line, (b) AG{x) diagram at various temper¬ 
atures. (c) Binodal and spinodal curves from TM (+ and -f), 
and MC simulations (*), compared with experimental data^® 
(— • ••) and the previous theoretical model (v) from Boukhris 
et al.®® 


B. The PbSi_i,Tea; alloy 

The atomic radius of Te is 24% larger than the S radius. 
Thus, the PbS-PbTe system follows the same trend as 
PbSe-PbTe, and they are immiscible at 0 K. 

Our AH values are greater than 0 eV for the whole 
range of concentrations (see Figure. 2(a)). CE predicted 



PbS JCjg PbTe 

FIG. 2. (a) Formation enthalpies of the PbS-PbTe structures 
using DFT calculations (O) a^nd CE technique (□). Highly 
symmetric structures are represented by A and the fitting of 
these points to obtain the interaction parameter is plotted with 
a blue dashed line, (b) AG(a;) diagram at various tempera¬ 
tures. (c) Binodal and spinodal curves from TM (+ and -f), 
and MC simulations (*), compared with experimental data®® 
(— ■ ••) and the previous theoretical model (v) from Boukhris 
et al.®® 


energies are again in agreement with DFT calculations, 
obtaining a CV score of about 3T0“^. The fitting for L{x) 
is also depicted in Figure. 2(a) using highly symmetric 
points. Similarly to PbSe-PbTe, the AG function changes 
to a convex shape at high temperatures (Figure. 2(b)). 

The calculated phase diagram of the PbSi_ 2 ,Tea; alloy 
is shown in Figure. 2(c). Experimental results show again 
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a slight asymmetry with a maximum around x = 0.3. 
This is in agreement with our results, while MC simu¬ 
lations fail to show this asymmetry. The predicted con- 
solute temperatures for PbS-PbTe follow the same trend 
as PbSe-PbTe. Results published by Boukhris et al.^^ 
considerably underestimate the experimental value for 
(1083 K). Our prediction of the consolute temperature is 
slightly larger than experiments^®; being 1385 and 1365 K 
using MC simulations and TM, respectively. As seen from 
Figure. 1(c) and Figure. 2(c) PbSe-PbTe and PbS-PbTe 
systems show a very similar trend of a slightly asymmet¬ 
ric spinodal curve, and considerably asymmetric binodal 
curve. This trend shows that formation of the Te-rich 
alloy starts at lower temperatures than Se-rich composi¬ 
tions. 


C. The PbSi-jjSej; alloy 

Our methodology was also applied to the PbS-PbSe 
system. All positive energies in Figure. 3 indicate that 
PbS and PbSe systems are not miscible at 0 K. As far as 
we know, there are no experimental data available below 
573 for this system. However, it has been shown that 
MC simulations predict the for different systems quite 
well, and can describe the miscibility gap.^^’®® For this al¬ 
loy, thermodynamic modelling predicts a Tc slightly below 
200 K and MC predicts a Tc slightly lower than 250 K. 


D. General considerations 


The lattice mismatch between the two solids (PbS, 
PbSe, or PbTe) is a good descriptor to analyze the trends 
observed experimentally (see Figure. 4). Lattice mis¬ 
match, e, is defined as: 


Solvent) 2^QQ 
^solvent 


( 11 ) 



PbS vsj. PbSe 


where &i^A,B)aCa denotes the lattice constant of intermedi¬ 
ate alloys and asoivent is the lattice constant of the most 
abundant binary alloy. There is a correlation between 
the lattice mismatch of the alloy, AiL, and the consolute 
temperature. The higher the mismatch, the higher A77 
becomes and thus, a higher Tc is obtained. For instance, 
the larger mismatch corresponds to the PbS-PbTe sys¬ 
tem, which presents a maximum enthalpy of formation 
at x = 0.5, with 80 meV/atom and a Tc of 1083 K. For 
the PbSe-PbTe system, the maximum AH is around 22 
meV/atom and the consolute temperature is 623 K. Fol¬ 
lowing this trend, PbS-PbSe system presents a smaller 
mismatch and a smaller AH 8 meV/atom. Thus, a con¬ 
solute temperature smaller than 623 K is expected. If 
we approximate this correlation to a linear function, for a 
lattice mismatch around 3% we get a Tc close to 270 K; 
which is in agreement with MC and our thermodynamic 
model results. 


FIG. 3. (a) Formation enthalpies of the PbS-PbSe structures 

using DFT calculations (O) and CE technique (□). Highly 
symmetric structures are represented by A and the fitting of 
these points to obtain the interaction parameter is plotted with 
a blue dashed line, (b) AG(a;) diagram at various tempera¬ 
tures. (c) Binodal and spinodal curves from TM (+ and -f), 
and MC simulations (*), compared with experimental data^® 
(— ■ •■) and the previous theoretical model (v) from Boukhris 
et al.^® 


Mismatch between lattices can be also used to explain 
the asymmetry of the binodal curves. We can define the 
asymmetry of the curve as the ratio between the decom¬ 
position temperature of two points equidistant to x = 0.5. 
We have chosen 0.2 and 0.8 to define our asymmetry de- 
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X 

FIG. 4. Lattice mismatch for lead chalcogenides alloys. For 
each ternary system we consider 0 mismatch when x = 0. 


scriptor, ct' 


T{x = 0.8) 

r(a; = 0.2)' 


( 12 ) 


Using this definition, we can assume that a perfectly sym¬ 
metric spinodal curve has eT=l. 

As discussed above, mismatch between lattices is di¬ 
rectly related to the magnitude and size of the spinodal 
curve. Mismatch is the driving force in the three systems 
we are studying. Thus, we propose a second asymmetry 
descriptor, e^., based on the ratio between the lattice mis¬ 
match at two points equidistant to a: = 0.5: 


IV. CONCLUSIONS 


A hybrid approach, comprising high-throughput ab- 
initio and cluster-expansion techniques is used to create 
a thermodynamic model for calculating binodal and spin¬ 
odal decompositions in pseudo binary lead chalcogenides 
(PbSe-PbTe and PbS-PbTe). The model overcomes the 
limitations of previous theoretical studies, where too few 
compositions and/or configurations were taken into ac¬ 
count. The obtained thermodynamical features are very 
close to the experimentally data, when available. We 
also capture the asymmetry of the binodal curve, exper¬ 
imentally observed and previously computationally unre¬ 
solved. Additionally, phase diagrams for systems without 
experimental characterization, such as the PbS-PbSe al¬ 
loy, are suggested. The results have been valitaded by 
using MC simulations, and lattice mismatch between the 
binary solids descriptors. 


Overall the approach is suitable for the high-throughput 
characterization of miscibility gaps, spinodal and other 
decomposition phenomena. 


e{x = 0.8) 
eix = 0 . 2 )' 


(13) 


The asymmetry descriptor values for the three systems 
are shown in Table I. er shows PbS-PbSe as the most 
asymmetric spinodal curve, then PbS-PbTe, and finally 
the PbSe-PbTe system. This trend is exactly the same 
for Cr, emphasizing the importance of the lattice strain in 
the of these systems. 


PbSe-PbTe PbS-PbTe PbS-PbSe 
(L86 080 080 

e,. 0.84 0.59 0.49 

TABLE 1. Asymmetric descriptors for spinodal curves for lead 
chacogenides 
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